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The Pairwise Einstein Model (PEM) of steps not only justifies the use of the Generalized Wigner 
Distribution (GWD) for Terrace Width Distributions (TWDs), it also predicts a specific form for 
the Step Position Distribution (SPD), i.e., the probability density function for the fluctuations of a 
step about its average position. The predicted form of the SPD is well approximated by a Gaussian 
with a finite variance. However, the variance of the SPD measured from either real surfaces or 
Monte Carlo simulations depends on Ay, the length of step over which it is calculated, with the 
measured variance diverging in the limit Ay oo. As a result, a length scale Lw can be defined 
as the value of Ay at which the measured and theoretical SPDs agree. Monte Carlo simulations of 
the terrace-step-kink model indicate that Lw ~ 14.2^q, where is the correlation length in the 
direction parallel to the steps, independent of the strength of the step-step repulsion. Lw can also 
be understood as the length over which a single terrace must be sampled for the TWD to bear a 
"reasonable" resemblence to the GWD. 

PACS numbers: PACS Number(s): 05.70.Np, 68.55.Jk, 68.35. Ct, 68.35.-p 



I. INTRODUCTION 

A key factor determining the equilibrium morphology 
of a vicinal crystal surface is the interaction between the 
steps on that surface. In many cases, the elastic and 
electronic contributions to the step-step interaction take 
the form 



V{L) 



A 

Z2 



(1) 



where A determines the strength of the step-step inter- 
action and L is the distance between steps. Because this 
is a typical step-step interaction, and because it has the 
remarkable property of yielding exact solutions to very 
plausible approximate theoriesi^, we confine ourselves 
in this paper to interactions of the form given in Eq. 
With this restriction, many of the quantities discussed in 
this paper depend only on a single dimensionless param- 
eter. 



A: 



f3A 



(2) 



where f3 is the step stiffness, /cb is Boltzmann's constant, 
and T is the absolute temperature. 

One of the easiest methods^i^ for experimentally de- 
termining the interaction between steps on a vicinal crys- 
tal surface is through the observation of the Terrace 
Width Distribution (TWD). Typically, this has been 
done by fitting the TWD to a Gaussian, which is a 
good approximation and justified by the Gruber-Mullins 
approximation2i& (analogous to the Einstein model^ of 
solids) if the steps strongly repel each other. The step- 
step interaction is then extracted from the variance of the 
Gaussian. Unfortunately, however, the Gaussian approx- 
imation is only good for strongly interacting steps, and 



there are conflicting theoriesSiSiiSiiiii^iiiii regarding the 
relationship between the step-step interaction and the 
variance. 

Over the past decade^i^ it has become apparent that 
the so-called Generalized Wigner Distribution (GWD) 
provides a much better approximation to the TWD. The 
GWD exhibits the positive skew observed in TWDs from 
experiments and simulations, and it is a good fit quan- 
titatively to TWDs produced from Monte Carlo simula- 
tions of the terrace-step-kink (TSK) model. More signifi- 
cantly, the GWD can be justified on the basis of plausible 



approximation 



.3.15 



the most important of which is that 



the interaction and fluctuations of two adjacent steps are 
explicitly considered; the Gruber-Mullins approximation 
only explicitly considers one step. The two steps are kept 
close to each other by a harmonic well, which approxi- 
mates the interactions with all other steps. This model is 
refered to as the Pairwise Einstein Model (PEM). Both 
the Gruber-Mullins and pairwise Einstein models start by 
interpreting the steps as world-lines of spinless fermions, 
with the y-direcion (the average direction of the steps) 
corresponding to time. 

This paper considers a different statistical measure of 
the vicinal surface: the Step Position Distribution (SPD). 
In Sec. ^ we show that the pairwise Einstein model 
predicts a Gaussian-like distribution for the position of 
steps. In Sec. lIIII these predictions are shown to compare 
well with numerical results from simulations of the TSK 
model, at least for systems of the "right size" . The depen- 
dence of the SPD on the length of the steps is discussed in 
Sec. IIVI for the purpose of comparison, the dependence 
of the TWD on the length of steps is likewise discussed 
in Sec. Finally, in Sec. I VII we summarize and draw 
our conclusions. 
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II. PREDICTIONS FROM THE PAIRWISE 
EINSTEIN MODEL 

As was shown in Ref. 0, the GeneraUzed Wigner Dis- 
tribution can be derived from a phenomenological treat- 
ment in which only two steps are treated expHcitly, the 
rest contributing a "confinement potential" related to the 
two-dimensional pressure and compressibility of the sys- 
tem of steps. We use the usual trick of mapping steps 
onto the worldlines of one-dimensional spinless fermions, 
which in this case have the Hamiltonianii^ 



In this dimensionless formulation, we require that 

{x2 - xi) = 1 ; 
this fixes the value of uj to 

oj — 2b„ 



where 



and 



= 1 + VI + 4I. 



After a change of variables^ii^ to 



Xc 



Xi + X2 



2 

S — X2 — Xi , 

this Hamiltonian becomes separabl o'^i^^ , 



(3) 
(4) 

(5) 
(6) 

(7) 

(8) 
(9) 



n 



and it has the remarkable property that all of the eigen- 
states are known. The only eigenstate of interest to us 
at present, however, is the ground state, which can be 
written^ii^ 



*o,o(s,a:cm) 



where 



2h 



: exp I 



(e+i)/2 



-'ibeX^^) 



r[(e + i)/2] 



(11) 



(12) 



is a constant of normalization. The probability of finding 
the combination a specific combination of relative separa- 
tion and "center of mass" is, of course, just q(s, Xcm), 



which can be rewritten in terms of the original variables 
xi and X2- 

P{xi,X2) = ^l^Q{s,Xcm) 

= -^{x2-xireM-2b,{xl+xl)U3) 



subject to the constraint X2>xi. We can integrate out 
all possible values of X2 to find the probability density 
function for xi: 



/>oo 

Qiixi) = / P{Xi,X2)dx2 



: exp(— 26^Xi) 



/•OO 

X / (x2 - xi)^exp(-26ga;2)dx2 . (14) 

As should be expected, the mean value of xi is —1/2 and 
the mean value of X2 is -M/2, so we define the analytic 
SPD to be the calculated probability density function for 
xi - {xi): 



Q(x) 



Qi ( a; + i 



: exp 



-2bn 
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1 



X2-x + -\ exp(-26ga;2) da;2(15) 



Although Q{x) can only be evaluated numerically (it 
can be rewritten as a complicated expression involving 
hypergeometric functions, but this does not seem to be 
genuinely helpful), it is straightforward, though tedious, 
to calculate its moments. The two most important are 
the mean, which is zero by definition, and the variance, 
which is given by 



Q 



2b„ 



1 



(16) 
(17) 



These two moments would be enough to entirely spec- 
ify the SPD if it were a Gaussian distribution, which it 
should be approximately; the Gruber-Mullins approxi- 
mation for the TWD, since it concerns the fiuctuations 
in position of only a single step, can be equally well inter- 
preted as an approximation for the SPD. In fact, both the 
coefficient of skewnesai^ and the kurtosisi^ of the theoret- 
ical SPD vanish in the limit of strong step-step repulsion. 
The coefficient of skewness is given asymptotically by 



71 



((xi - {x^)f) V6 , 



/2. 



18 



(18) 



note that the coefficient of skewness would have the op- 
posite sign if it had been defined as ((x2 — {x2)y) crg^-y^ . 
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The kurtosis, which is the same regardless of which step 
is considered, is given asymptoticaUy by 



72 



((^1 - (^i»"> 



1 

12' 



(19) 



The fact that the kurtosis is not exactly zero is not in 
itself surprising; even within the Gruber-MuUins approx- 
imation, the Gaussian distribution is only obtained in 
the limit of large A. The symmetry of our original prob- 
lem of an infinite number of steps on an infinite vicinal 
surface, on the other hand, means that the coefficient 
of skewness, by contrast, must be zero for the original 
problem. Any given step on the surface can be consid- 
ered "step 1" , with its downhill neighbor as "step 2" , or 
it can be considered "step 2" , with its uphill neighbor as 
"step 1" ; calling it one or the other breaks the symmetry 
and permits a nonzero coefficient of skewness. 
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FIG. 1: An illustration of problem that can be caused by 
the variability of the mean step position when the step posi- 
tion distribution (SPD) is calculted from numerical or exper- 
imental results. In this example, Gaussian distributions with 
identical variances ((t^ = 2.5^) are binned into histograms by 
means of Eq. 1201 . The only differences between the three 
distributions are the values of fj,: circles, ^ — 0.452; squares, 
/i=— 0.279; diamonds, fi — —0.131. 



III. COMPARISON WITH MONTE CARLO 
SIMULATIONS 

In order to test the applicability of Eq. (|15|l . we have 
performed Monte Carlo simulations of the terrace-step- 
kink (TSK) model and measured the SPD for several 
values of A. 

The geometry of the simulated systems was as follows. 
All simulations were for systems of 30 steps; the length 
of each of which was Ly = 1000a (where a is the lat- 
tice constant) in the average direction of the steps (the 
y-direction in "Maryland notation" ) . The mean step sep- 
aration was (L) = 10a, and periodic boundary conditions 
were applied. 

The dynamic used was a local Metropolis update. The 
temperature was set at = O.Abe, where e is the kink 
energy; in a previous study, this was approximately the 
temperature at which TWDs from the restricted TSK 
model showed the best agreement with the Generalized 
Wigner Distribution. Each simulation was equilibrated 
for at least 500 000 Monte Carlo steps per site (MCSS) at 
the temperature and value of A at which measurements 
were taken; the initial configurations, however, were not 
typically straight steps, but steps that had been equili- 
brated at some other value of A. Data were taken from 
1 000 "snapshots," taken at intervals of 1 000 MCSS. 

Although the terrace width is always an integer multi- 
ple of a in the TSK model, the average step position can 
be any rational number, depending only on the size of 
the simulation. Since the step position x is always an in- 
teger, the histogram of positions for any given step need 
not be symmetric. 

In order to show concretely what this means, consider 
a situation in which a Gaussian distribution with mean 
/i and variance is binned into a histogram as follows. 
The weight assigned to each integer k is given by inte- 



grating the Gaussian between fc— 1/2 and fc-t-1/2: 



W{k) 



erf 



erf 



fc+l/2 



exp 



fc-l/2 



(x-M)' 
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dec 



2a 

fc + (1/2) 



2ct 



(20) 



For our example, we choose a = 2.5 and three "random" 
values of n between -0.5 and +0.5. The results are shown 
in Fig. n Clearly none of the histograms is completely 
symmetric, and the differences between them are note- 
worthy. 

Something similar can and does happen when the 
SPDs are calculated from Monte Carlo simulations by 
binning the positions into histograms. As a result, the 
statistical uncertainties are considerably larger than they 
are for the corresponding TWDs, and the SPDs are not 
perfectly symmetric about their peaks, as can be seen in 
Figs. 121 and 13 Note the qualitative similarities between 
the Monte Carlo results (circles) in Figs. |21 and |21 and the 
values of W{k) for /i = — 0.279 (squares) and /x = — 0.131 
(diamonds) in Fig.Q] This agreement suggests that dur- 
ing the process of equlibration, the majority of the steps 
moved slightly to the left {i.e., uphill). 

In spite of this, the agreement of the SPDs calcu- 
lated from simulations and the theoretical Q{x) calcu- 
lated from Eq. (|15|) is reasonably good. Even more im- 
pressive is the agreement between Q{x) and the Gaussian 
with zero mean and variance given by Eq. (|16|l . Although 
Eqs. H18() and H19() suggest that the Gaussian approxima- 
tion will be increasingly good as A becomes large, it is 
clear from the figures that the Gaussian approximation 
is good for even for A = 0. 
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FIG. 2: Comparison of the SPD for i = given by Eq. lO 
(solid curve) with a histogram SPD from a Monte Carlo sim- 
ulation (symbols). Also shown is a Gaussian (dotted curve) 
with a me; '' ' ' • i y ,.j 




FIG. 3: Comparison of the SPD for i = 8 given by Eq. ^ 
(solid curve) with a histogram SPD from a Monte Carlo sim- 
ulation (symbols). Also shown is a Gaussian (dotted curve) 
with a mean of zero and a variance given by Eq. 1161 . 



IV. SCALING OF THE SPD 

Although the agreement between Eq. (|15|l and the nu- 
merical SPDs discussed above is highly suggestive, it is 
clear that actual step position distributions must depend 
on the length Ay of step over which they are averaged. 
This is best demonstrated by considering the variance of 
a measured SPD, which is given by 



4(Ay) ^ {Ay)-' ( / [x{y) 

\ J-Ay/2 



where 



(Ay)- 



Ay/2 



Ay/2 



x{y)dy. 



(21) 



(22) 



(23) 



which characterizes the wandering of an individual 

Stepl7,18,19,20,21_ 



Clearly, aQ{Ay) is closely related t( 



g.,{Ay) = ([x{Ay) - x{0)Y 



FIG. 4: Comparison of the correlation length calculated from 
the SPD (^q) and from the Gruber-MuUins approximation 
(Cgm), evaluated numerically for (L) — 10 and ksT = 0.45e. 
Although there is decent agreement for A> 5, ^gm unphys- 
ically diverges as j4 — > 0. In contrast, remains finite and 
reasonable for all nonegative values of A. 



It is tempting to identify x, the average value of x for 
a particular conformation of a step, with x(0), the value 
of X at the average y-position. This leads to 



a^QiAy) « (Ay)- 
For small Ay, « ci|Ay^ 



Ay/2 



Ay/2 



9x{y)dy. 



(24) 



17.18.19.20.21. 



Eq. if^ implies 



i(ci/2)Ay. Likewise, for large Ay,i7ii8,i9,20,2i 
9x{Ay) « C2 + czln\Ay\ 



and Eq. implies 



)(Ay) 



C4 -I- c^ln\Ay\ 



(25) 



(26) 



The observation, made in the previous section, that 
Q{x) is to a very good approximation Gaussian is help- 
ful towards the calculation of the characteristic length 
for cr|. In Ref.H the "TWD" was calculated within the 
Gruber-MuUins approximation; because the position of 
only one step was explicitly taken into account, though, 
it could equally be considered a SPD. In fact, the Gaus- 
sian solution is a more appropriate description of a SPD, 
which is symmetric, than a TWD, which is asymmetric. 
Substituting the variance of the SPD into Eq. (18) of 
Ref. 1^ we find the correlation length to be 



,w 



fcRT 



(27) 



Figure 0| shows a comparison between and the corre- 
lation length from Ref. 0. 

Scaled by ctq ^'^^ ^Q' '^^{^u) appears to be inde- 
pendent of A\ although the PEM incorrectly predicts that 
crQ(Ay) remains finite in the limit Ay— >oo, it neverthe- 
less provides the correct scaling factors. This is remark- 
able, since although ^^.(Ay) shows scaling with temper- 
ature, it does not exhibit scaling independent of 
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FIG. 5: A power-law fit to all the Monte Carlo estimates of 
(7q(Aj/) for < A < 8, Ay < ^q, indicates an initial growth 
with an exponent of 0.797 ± 0.006. 




FIG. 7: Terrace Width Distributions calculated between a 
single pair of neighboring steps depend on Ay, the length of 
step over which the distribution is averaged. Although in the 
limit Ay — > oo the TWD converges to the Generalized Wigner 
Distribution (to a very good approximation), when Ay/^Q is 
small the TWD is dominated by noise. These results are 
typical for ^ = 0. 



find a "Wigner length", Lw, defined by 



to be 



(14.183 ±0.035)^Q. 



(28) 



(29) 



SCALING OF THE TWD 



FIG. 6: A fit to the Monte Carlo estimates of a'Q(Ay) for 
A = 8, 4^(3 <Ay<Ly/2, indicates an asymptotic growth given 
by o-|(At/)/(7^,w~ 0.158+0.3175 ln(Ay/^Q). The length Lw, 
defined by Eq. o-q(I/w) = cr^, is consequently given 

by Lw = (14.183 ±0.035)Cq. 



For Ay < ^q, a least-squares fit indicates power-law 
growth of aqlAy) with an exponent of 0.797 ± 0.006 
(Fig. O . Equation H24|l predicts power-law growth, but 
with an exponent of 1. Interestingly, the power-law be- 
havior of gx{Ay) extends only out toi^' Ay « O.lj/coii; 
since j/coii = Cq/I'"' ~ 2) (for A = 0), power-law scaling 
extends farther for (Tg(Ay) than for ga:{Ay). 

For large Ay, aQ{Ay) follows the logarithmic scaling 
of Eq. 1)26(1 . A least-squares fit was performed on the 
A = 8 data, since this has the smallest value of among 
the available simulations, and hence the largest avail- 
able values of Ay/^g. To avoid the crossover from the 
power-law regime, the fit was restricted to Ay > 4^q; 
likewise, the fit was limited to Ay<Ly/2 to limit finite- 
size effects. The resulting fit, shown in Fig. is in good 
agreement with data from all values of A except where 
finite-size effects become evident. The fitted parameters, 
C4 = 0.1578 ± 0.0004 and C5 = 0.3175 ± 0.0002, allow us to 



It seems somewhat surprising that so many correlation 
lengths are necessary for the PEM to agree with the ob- 
served variance. To better understand this, it is helpful 
to consider the corresponding scaling of the TWD when 
it is calculated under the same restrictions as crQ(Ay). 
Specifically, the TWD must be averaged over a given 
length. Ay, of a single pair of adjacent steps in a sin- 
gle "snapshot". This is very different from the analysis 
presented in Ref. l22l where, as in other previous work, 
averages were made over the entire length Ly of the sim- 
ulations, over all pairs of neighboring steps, and over all 
"snapshots". Remembering that the j/-direction corre- 
sponds to time in the worldline interpretation of steps, 
the averages we are about to calculate correspond to 
time averages in statistical mechanics, whereas the pre- 
vious averages have combined the time average with two 
kinds of ensemble average (over different pairs and dif- 
ferent "snapshots"). Only in the limits of long times and 
large ensembles should one expect these averages to be 
identicals^. 

In the language of Ref. 17, £,q is approximately the dis- 
tance between "collisions" of neighboring steps. In order 
to sample the distribution of terrace widths adequately, 
a step must "collide" several times with its neighbors. 
This is shown clearly in Figs. and |H1 For the simula- 
tion parameters given in Sec. IIIII Lw > Ly for A — 0; 
consequently, the TWDs shown in Fig. [7| are dominated 
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FIG. 8: For small Ay/£,Q, typical TWDs are dominated by 
noise, but for Ay > Lw ~ 14.25q, the Generalized Wigner 
Distribur ^ • ■ - . x o 




15 20 25 30 



FIG. 9: The variance of the TWD also approximately scales 
with Ay/^Q. At Ay = Lw ~ 14.24q, the average variance 
of TWDs generated from single pairs of neighboring steps is 
within about %10 of the variance given by the GWD. 



by noise. For A = 8, on the other hand, Lw < Ly, and 
we are able to see in Fig.|Slthe crossover into the regime 
where the GWD, not noise, is dominant. 

The role of step collisions in equilibrating the TWD can 
also be seen from ap{Ay), which is the variance of TWDs 
calculated from a length Ay of neighboring steps, aver- 
aged over all pairs of neighboring steps, starting points, 
and "snapshots". Like aQ{Ay), hmAy^o CTp(Aj/) — 0; 
unlike (TQ(Ay), limAj,-»oo c''p(A?/) is finite and given 
approximateljiS^ by the PEM resulli^. 



'p.w 



g+l 
'2b7 



(30) 



This suggests plotting 1 — tTp(Aj/)/(Tp^ vs. Ay/£_Q to 
determine whether the approach to the asymptotic limit 
is exponential or power-law. As shown in Fig. El the 
scaling appear to be neither a simple power law nor a 
simple exponential decay, but it is difficult to be certain 
since the TWD does not converge exactly to the GWD 
even in the limit Ay oo. Also, the scaling does not 
appear to be quite as precise as in Figs. |31andEl this 
is not surprising, since the correlation length ^p for the 



TWD is not identical to ^q. 

More significantly. Fig. El indicates that crp(Lw) is 
within about %10 of the approximate asymptotic value, 
CTp^. This is a very plausible threshold for statistics 
from the TWD. 



VI. CONCLUSION 

For the common case in which steps on a vicinal crys- 
tal surface interact according to Eq. (Q, the generalized 
Wigner distribution (GWD) has been shown previously 
to be in excellent agreement with the terrace width dis- 
tribution (TWD). To fully appreciate the model which 
predicts the GWD, though, it is necessary to examine its 
predictions for other statistical properties and how well 
these predictions agree with actual measurements. This 
article has made such a comparison between the pre- 
dicted and measured step position distribution (SPD). 
The results demonstrate both the strength and limita- 
tions of the Pairwise Einstein Model (PEM). 

Since the SPD is so well approximated by a Gaussian, 
it is tempting to compare it directly with Gaussian theo- 
ries of the TWD. As can be seen in Table 1, in the limit 
of strongly interacting steps the variance of the SPD is 
slightly larger than that of the Gruber-MuUins approx- 
imation, but less than the variance of the TWD given 
by either the "Saclay" or "modified Grenoble" approxi- 
mations. This is reasonable; unlike the Gruber-MuUins 
Hamiltonian, Eq. Q does not have fixed walls, so the 
steps can experience larger fluctuations. In spite of this, 
since the Gruber-MuUins approximation allows only one 
step to move, it can be regarded equally as an approxima- 
tion for the TWD or for the SPD. The fact that the SPD 
is smaller than the other approximations of the TWD 
is apparently due to the fact that correlations between 
fluctuations of adjacent steps are to some degree taken 
into account in all these approximations, so that they are 
specifically approximations for the TWD, not the SPD. 



Because the PEM confines both steps within a har- 
monic well, the theoretical asymptotic variance of the 
SPD must be finite. However, the vicinal surface is rough, 
and the variance of the SPD diverges logarithmically with 
the length of step Ay from which it is calculated. At some 
finite length, Lw, the prediction of the PEM is accurate. 
As Fig.jSlshows, Lw~14.2^q. That so many "collisions" 
between neighboring steps are needed to adequately sam- 
ple the statistics resulting from their interactions is sup- 
ported by obeservations of the dependence of the TWD 
on Ay, as shown in Figs.[71|Hl and El 

In principle, the SPD could be used to determine A. 
However, the SPD is strongly affected by the random 
position of the average step position, and it depends far 
too strongly on Ay. The TWD has neither of these two 
restrictions and is a more practical alternative for deter- 
mining A. Instead, the utility of the SPD lies in clarifying 
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TABLE I: Asymptotic variances in the limit of strong step- 
step repulsion. The Gaussian-like approximation for the step 
position distribution (SPD) given by Eq. 1151 is compared 
with selected approximations for the TWD. Except for the 
Generalized Wigner Distribution, all approximate TWDs are 
Gaussian approximations. Note also that our approximation 
for the SPD and the Generalized Wigner Distribution are 
both independent of the number of interacting steps, whereas 
the Gaussian approximations are not. (See also Table 1 of 
Ref^.) 

Distribution Reference Asymptotic Variance 



SPD Eq. JUJ 0.375q-^ 

Generalized Wigner 24^2^ 0.5^"^ 

Gruber Mullins (aU steps) 0.278^"^ 

" (nearest neighbors) " 0.289g~^ 

Modified Grenoble (aU steps) iSJi^ 0.495£|-^ 

" (nearest neighbors) " 0.520£i~ 

Saclay (all steps) 0.405g~^ 



the Pairwise Einstein Model. The finite length L-yv in- 
troduced in this work emerges more naturally than the 
finite length V that was introduced in Refs. and 0, 
but the two are obviously related. Both help describe a 
short-lived dynamic constraint that is roughly analagous 
to a reptation tubc^^ in polymer physics. 

Naturally, the remarkable success of the Pairwise Ein- 
stein Model suggests that a Debye model^ might lead to 
even better descriptions of vincinal crystal surfaces. Pre- 
liminary results^ from such studies correctly show that 
gx{Ay) diverges logarithmically. 
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